Molecular Simulation of Argon Adsorption and Diffusion in a Microporous Carbon with Poroelastic Couplings.

Neglected for a long time in molecular simulations of fluid adsorption and transport in microporous carbons, adsorption-induced deformations of the matrix have recently been shown to have important effects on both sorption isotherms and diffusion coefficients. Here we investigate in detail the behavior of a recently proposed 3D-connected mature kerogen model, as a generic model of aromatic microporous carbon with atomic H/C ∼ 0.5, in both chemical and mechanical equilibrium with argon at 243 K over an extended pressure range. We show that under these conditions the material exhibits some viscoelasticity, and simulations of hundreds of nanoseconds are required to accurately determine the equilibrium volumes and sorption loadings. We also show that neglecting matrix internal deformations and swelling can lead to underestimations of the loading by up to 19% (swelling only) and 28% (swelling and internal deformations). The volume of the matrix is shown to increase up to about 8% at the largest pressure considered (210 MPa), which induces an increase of about 33% of both pore volume and specific surface area via the creation of additional pores, yet does not significantly change the normalized pore size distribution. Volume swelling is also rationalized by using a well-known linearized microporomechanical model. Finally, we show that self-diffusivity decreases with applied pressure, following an almost perfectly linear evolution with the free volume. Quantitatively, neglecting swelling and internal deformations tends to reduce the computed self-diffusivities.

Texture, Nanotexture, and Structure of Carbon Nanotube-Supported Carbon Cones.

Graphene-based carbon micro-/nano-cones were prepared by depositing pyrolytic carbon onto individual carbon nanotubes as supports using a specific chemical vapor deposition process. They were investigated by means of high-resolution scanning electron microscopy, low-voltage aberration-corrected transmission electron microscopy, Raman spectroscopy, and molecular dynamics modeling. While the graphenes were confirmed to be perfect, the cone texture was determined to be preferably scroll-like, with the scroll turns being parallel to the cone axis. Correspondingly, many of the concentrically displayed graphenes (actually scroll turns) exhibit the same helicity vector. When radii of curvature are large enough, this could allow for coherent stacking to locally take place in spite of the lattice shift induced by the curvature. A particular care was taken on investigating the cone apexes, in which a specific type of graphene termination was observed, here designated as the "zip" defect. Calculations determined a plausible stable structure that such a defect type may correspond to. This defect was found to generate a very low Raman ID/ID' band ratio (1.5), for which physical reasons are proposed. Combining our results and that of the literature allowed proposing an identification chart for a variety of defects able to affect the graphene lattice or edges.

Timescale prediction of complex multi-barrier pathways using flux sampling molecular dynamics and 1D kinetic integration: Application to cellulose dehydration.

Reactive molecular dynamics (MD) simulations, especially those employing acceleration techniques, can provide useful insights on the mechanism underlying the transformation of buried organic matter, yet, so far, it remains extremely difficult to predict the time scales associated with these processes at moderate temperatures (i.e., when such time scales are considerably larger than those accessible to MD). We propose here an accelerated method based on flux sampling and kinetic integration along a 1D order parameter that can considerably extend the accessible time scales. We demonstrate the utility of this technique in an application to the dehydration of crystalline cellulose at temperatures ranging from 1900 K to 1500 K. The full decomposition is obtained at all temperatures apart from T = 1500 K, showing the same distribution of the main volatiles (H2O, CO, and CO2) as recently obtained using replica exchange molecular dynamics. The kinetics of the process is well fitted with an Arrhenius law with Ea = 93 kcal/mol and k0 = 9 × 1019 s-1, which are somehow larger than experimental reports. Unexpectedly, the process seems to considerably slow down at lower temperatures, severely departing from the Arrhenius regime, probably because of an inadequate choice of the order parameter. Nevertheless, we show that the proposed method allows considerable time sampling at low temperatures compared to conventional MD.

Methane Diffusion in a Flexible Kerogen Matrix.

It has been recognized that the microporosity of shale organic matter, especially that of kerogen, strongly affects the hydrocarbon recovery process from unconventional reservoirs. So far, the numerical studies on hydrocarbon transport through the microporous phase of kerogen have neglected the effect of poromechanics, that is, the adsorption-induced deformations, by considering kerogen as a frozen, nondeformable, matrix. Here, we use molecular dynamics simulations to investigate methane diffusion in an immature (i.e., with high H/C ratio) kerogen matrix, while explicitly accounting for adsorption-induced swelling and internal matricial motions, covering both phonons and nonperiodic internal deformations. However, in the usual frozen matrix approximation, diffusivity decreases with increasing fluid loading, as evidenced by a loss of free volume, accounting for adsorption-induced swelling that gives rise to an increase in free volume and, hence, in diffusivity. The obtained trend is further rationalized using a Fujita-Kishimoto free volume theory initially developed in the context of diffusion in swelling polymers. We also quantify the enhancing effect of the matrix internal motions (i.e., at fixed volume) and show that it roughly gives an order of magnitude increase in diffusivity with respect to a frozen matrix, thanks to fluctuations in the pore connectivity. We eventually discuss the possible implications of this work to explain the productivity slowdown of hydrocarbon recovery from shale immature reservoirs.

Poroelasticity of Methane-Loaded Mature and Immature Kerogen from Molecular Simulations.

While hydrocarbon expulsion from kerogen is certainly the key step in shale oil/gas recovery, the poromechanical couplings governing this desorption process, taking place under a significant pressure gradient, are still poorly understood. Especially, most molecular simulation investigations of hydrocarbon adsorption and transport in kerogen have so far been performed under the rigid matrix approximation, implying that the pore space is independent of pressure, temperature, and fluid loading, or in other words, neglecting poromechanics. Here, using two hydrogenated porous carbon models as proxies for immature and overmature kerogen, that is, highly aliphatic hydrogen-rich vs highly aromatic hydrogen-poor models, we perform an extensive molecular-dynamics-based investigation of the evolution of the poroelastic properties of those matrices with respect to temperature, external pressure, and methane loading as a prototype alkane molecule. The rigid matrix approximation is shown to hold reasonably well for overmature kerogen even though accounting for flexibility has allowed us to observe the well-known small volume contraction at low fluid loading and temperature. Our results demonstrate that immature kerogen is highly deformable. Within the ranges of conditions considered in this work, its density can double and its accessible porosity (to a methane molecule) can increase from 0 to ∼30%. We also show that these deformations are significantly nonaffine (i.e., nonhomogeneous), especially upon fluid adsorption or desorption.

From cellulose to kerogen: molecular simulation of a geological process.

The process by which organic matter decomposes deep underground to form petroleum and its underlying kerogen matrix has so far remained a no man's land to theoreticians, largely because of the geological (Myears) timescale associated with the process. Using reactive molecular dynamics and an accelerated simulation framework, the replica exchange molecular dynamics method, we simulate the full transformation of cellulose into kerogen and its associated fluid phase under prevailing geological conditions. We observe in sequence the fragmentation of the cellulose crystal and production of water, the development of an unsaturated aliphatic macromolecular phase and its aromatization. The composition of the solid residue along the maturation pathway strictly follows what is observed for natural type III kerogen and for artificially matured samples under confined conditions. After expulsion of the fluid phase, the obtained microporous kerogen possesses the structure, texture, density, porosity and stiffness observed for mature type III kerogen and a microporous carbon obtained by saccharose pyrolysis at low temperature. As expected for this variety of precursor, the main resulting hydrocarbon is methane. The present work thus demonstrates that molecular simulations can now be used to assess, almost quantitatively, such complex chemical processes as petrogenesis in fossil reservoirs and, more generally, the possible conversion of any natural product into bio-sourced materials and/or fuel.

Mechanism of strength reduction along the graphenization pathway.

Even though polycrystalline graphene has shown a surprisingly high tensile strength, the influence of inherent grain boundaries on such property remains unclear. We study the fracture properties of a series of polycrystalline graphene models of increasing thermodynamic stability, as obtained from a long molecular dynamics simulation at an elevated temperature. All of the models show the typical and well-documented brittle fracture behavior of polycrystalline graphene; however, a clear decrease in all fracture properties is observed with increasing annealing time. The remarkably high fracture properties obtained for the most disordered (less annealed) structures arise from the formation of many nonpropagating prefracture cracks, significantly retarding failure. The stability of these reversible cracks is due to the nonlocal character of load transfer after a bond rupture in very disordered systems. It results in an insufficient strain level on neighboring bonds to promote fracture propagation. Although polycrystallinity seems to be an unavoidable feature of chemically synthesized graphenes, these results suggest that targeting highly disordered states might be a convenient way to obtain improved mechanical properties.

Rippled nanocarbons from periodic arrangements of reordered bivacancies in graphene or nanotubes.

We report on various nanocarbons formed from a unique structural pattern containing two pentagons, three hexagons, and two heptagons, resulting from local rearrangements around a divacancy in pristine graphene, or nanotubes. This defect can be inserted in sheets or tubes either individually or as extended defect lines. Sheets or tubes containing only this defect as a pattern can also be obtained. These fully defective sheets, and most of the tubes, present a very pronounced rippled (wavy) structure and their energies are lower than other structures based on pentagons and heptagons published so far. Another particularity of these rippled carbon sheets is their ability to fold themselves into a two-dimensional porous network of interconnected tubes upon heat treatment as shown by hybrid Monte Carlo simulations. Finally, contrary to the common belief that pentagon/heptagon based structures are metallic, this work shows that this defect pattern should give rise to semimetallic conduction.

Reaction mechanism for the thermal decomposition of BCl3/CH4/H2 gas mixtures.

This paper presents an ab initio study of the B/C/Cl/H gas phase mechanism, featuring 10 addition-elimination reactions involving BH(i)Cl(j) (i + j ≤ 3) species and a first description of the chemical interaction between the carbon-containing and boron-containing subsystems through the three reactions BCl(3) + CH(4) ⇌ BCl(2)CH(3) + HCl, BHCl(2) + CH(4) ⇌ BCl(2)CH(3) + H(2), and BCl(2) + CH(4) ⇌ BHCl(2) + CH(3). A reaction mechanism is then proposed and used to perform some illustrative equilibrium and kinetic calculations in the context of chemical vapor deposition (CVD) of boron carbide. Our results show that the new addition-elimination reaction paths play a crucial role by lowering considerably the activation barrier with respect to previous theoretical evaluations; they also confirm that BCl(2)CH(3) is an important species in the mechanism.

Theoretical study of the decomposition of BCl3 induced by a H radical.

We report on a theoretical study of the gas-phase decomposition of boron trichloride in the presence of hydrogen radicals using ab initio energetic calculations coupled to TST, RRKM, and VTST-VRC kinetic calculations. In particular, we present an addition-elimination mechanism (BCl(3) + H → BHCl(2) + Cl) allowing for a much more rapid consumption of BCl(3) than the direct abstraction reaction (BCl(3) + H → BCl(2) + HCl) considered up to now. At low temperatures, T ≤ 800 K, our results show that a weakly stabilized complex BHCl(3) is formed with a kinetic law compatible with the consumption rate measured in the former experiments. At higher temperatures, this complex is not stable and then easily eliminates a chlorine atom. Our work also shows that a very similar mechanism, involving the same intermediate and sharing the same transition state, allows for the elimination of HCl. A dividing coefficient between these two elimination pathways is obtained from both a potential energy surface based statistical analysis and an ab initio molecular dynamics transition path sampling simulation. It finally allows partitioning of the global consumption rate of BCl(3) in terms of the formation of (i) BHCl(3), (ii) BHCl(2) + Cl through a H addition/Cl elimination mechanism, (iii) BCl(2) + HCl through a H addition/HCl elimination mechanism, and (iv) BCl(2) + HCl through direct abstraction.

Hindered rotor models with variable kinetic functions for accurate thermodynamic and kinetic predictions.

We present an extension of some popular hindered rotor (HR) models, namely, the one-dimensional HR (1DHR) and the degenerated two-dimensional HR (d2DHR) models, allowing for a simple and accurate treatment of internal rotations. This extension, based on the use of a variable kinetic function in the Hamiltonian instead of a constant reduced moment of inertia, is extremely suitable in the case of rocking/wagging motions involved in dissociation or atom transfer reactions. The variable kinetic function is first introduced in the framework of a classical 1DHR model. Then, an effective temperature and potential dependent constant is proposed in the cases of quantum 1DHR and classical d2DHR models. These methods are finally applied to the atom transfer reaction SiCl(3)+BCl(3)→SiCl(4)+BCl(2). We show, for this particular case, that a proper accounting of internal rotations greatly improves the accuracy of thermodynamic and kinetic predictions. Moreover, our results confirm (i) that using a suitably defined kinetic function appears to be very adapted to such problems; (ii) that the separability assumption of independent rotations seems justified; and (iii) that a quantum mechanical treatment is not a substantial improvement with respect to a classical one.

Molecular dynamics of carbon dioxide, methane and their mixtures in a zeolite possessing two independent pore networks as revealed by computer simulations.

The molecular motion of methane (CH(4)) and carbon dioxide (CO(2)) sorbed in the two independent pore networks, being termed hereafter as large cavity (LC) and sinusoidal channel (SC) regions of the siliceous MWW-framework-type zeolite ITQ-1, is studied by means of atomistic computer simulation. Equilibrium molecular dynamics predicts different loading dependences for the two gases, for both the self and the collective (Maxwell-Stefan) diffusion coefficients; in particular, the transport coefficients of CH(4) go through a maximum as its loading in the zeolite increases, whereas CO(2) dynamics exhibits the decreasing trend with loading usually observed in nanoporous materials. The different loading dependence of the self-diffusivity for the two sorbates is attributed to their different probability density distribution within the supercages in the LC system of the ITQ-1 unit cell. The composition and occupancy dependence of the self-diffusivity of each component in their binary mixtures can be explained in terms of the selectivity for CO(2) sorption thermodynamics in the zeolite. The collective diffusivity loading dependence of the single and binary sorbate system is explainable on the basis of the strength of intermolecular interactions along the diffusion direction connecting the supercages by invoking the quasichemical mean field theory.

The self-referential method for linear rigid bodies: application to hard and Lennard-Jones dumbbells.

The self-referential (SR) method incorporating thermodynamic integration (TI) [Sweatman et al., J. Chem. Phys. 128, 064102 (2008)] is extended to treat systems of rigid linear bodies. The method is then applied to obtain the canonical ensemble Helmholtz free energy of the alpha-N(2) and plastic face centered cubic phases of systems of hard and Lennard-Jones dumbbells using Monte Carlo simulations. Generally good agreement with reference literature data is obtained, which indicates that the SR-TI method is potentially very general and robust.

The self-referential method combined with thermodynamic integration.

The self-referential method [M. B. Sweatman, Phys. Rev. E 72, 016711 (2005)] for calculating the free energy of crystalline solids via molecular simulation is combined with thermodynamic integration to produce a technique that is convenient and efficient. Results are presented for the chemical potential of hard sphere and Lennard-Jones face centered cubic crystals that agree well with this previous work. For the small system sizes studied, this technique is about 100 times more efficient than the parameter hopping technique used previously.

Hit and miss of classical nucleation theory as revealed by a molecular simulation study of crystal nucleation in supercooled sulfur hexafluoride.

Classical nucleation theory pictures the homogeneous nucleation of a crystal as the formation of a spherical crystalline embryo, possessing the properties of the macroscopic crystal, inside a parent supercooled liquid. In this work we study crystal nucleation in moderately supercooled sulfur hexafluoride by umbrella sampling simulations. The nucleation free energy evolves from 5.2kBT at T=170 K to 39.1kBT at T=195 K. The corresponding critical nucleus size ranges from 40 molecules at T=170 K to 266 molecules at T=195 K. Both nucleation free energy and critical nucleus size are shown to evolve with temperature according to the equations derived from the classical nucleation theory. Inspecting the obtained nuclei we show, however, that they present quite anisotropic shapes in opposition to the spherical assumption of the theory. Moreover, even though the critical nuclei possess the structure of the stable bcc plastic phase, the only mechanically stable crystal phase for SF6 in the temperature range investigated, they are shown to be less ordered than the corresponding macroscopic crystal. Their crystalline order is nevertheless shown to increase regularly with their size. This is confirmed by a study of a nucleus growth from a critical size to a size of the order of 10(4) molecules. Similarly to the fact that it does not affect the temperature dependence of the nucleation free energy and of the critical nucleus size, the ordering of the nucleus with size does not affect the growth rate of the nucleus.

Sorption thermodynamics of CO2, CH4, and their mixtures in the ITQ-1 zeolite as revealed by molecular simulations.

The thermodynamic properties and siting of carbon dioxide and methane sorbed in the siliceous form of zeolite MCM-22, ITQ-1, were studied by means of grand canonical Monte Carlo simulation. ITQ-1 comprises two independent pore systems of different geometry. It was found to be CO(2)-selective toward CO(2)/CH(4) gas mixtures, its equilibrium selectivity being distinctly higher in its sinusoidal channel pore system than in the large cavity system over a wide range of pressures starting from the Henry law regime, at the three temperatures considered. A maximum in selectivity is observed at low temperature, high pressure, and methane-rich gas-phase composition.

Molecular simulation of the homogeneous crystal nucleation of carbon dioxide.

We report on a molecular simulation study of the homogeneous nucleation of CO2 in the supercooled liquid at low pressure (P = 5 MPa) and for degrees of supercooling ranging from 32% to 60%. In all cases, regardless of the degree of supercooling, the structure of the crystal nuclei is that of the Pa3 phase, the thermodynamically stable phase. For the more moderate degree of supercooling of 32%, the nucleation is an activated process and requires a method to sample states of high free energy. In this work, we apply a series of bias potentials, which promote the ordering of the centers of mass of the molecules and allow us to gradually grow crystal nuclei. The reliability of the results so obtained is assessed by studying the evolution of the nuclei in the absence of any bias potential, and by determining their probability of growth. We estimate that the size of the critical nucleus, for which the probability of growth is 0.5, is approximately 240 molecules. Throughout the nucleation process, the crystal nuclei clearly exhibit a Pa3 structure, in apparent contradiction with Ostwald's rule of stages. The other polymorphs have a much larger free energy. This makes their formation highly unlikely and accounts for the fact that the nucleation of CO2 proceeds directly in the stable Pa3 structure.

Atomistic simulation of the homogeneous nucleation and of the growth of N2 crystallites.

We report on a computer simulation study of the early stages of the crystallization of molecular nitrogen. First, we study how homogeneous nucleation takes place in supercooled liquid N(2) for a moderate degree of supercooling. Using the umbrella sampling technique, we determine the free energy barrier of formation for a critical nucleus of N(2). We show that, in accord with Ostwald's rule of stages, the structure of the critical nucleus is predominantly that of a metastable polymorph (alpha-N(2) for the state point investigated). We then monitor the evolution of several critical nuclei through a series of unbiased molecular dynamics trajectories. The growth of N(2) crystallites is accompanied by a structural evolution toward the stable polymorph beta-N(2). The microscopic mechanism underlying this evolution qualitatively differs from that reported previously. We do not observe any dissolution or reorganization of the alpha-like core of the nucleus. On the contrary, we show that alpha-like and beta-like blocks coexist in postcritical nuclei. We relate the structural evolution to a greater adsorption rate of beta-like molecules on the surface and show that this transition actually starts well within the precritical regime. We also carefully investigate the effect of the system size on the height of the free energy barrier of nucleation and on the structure and size of the critical nucleus.

Reorganization and growth of metastable alpha-N2 critical nuclei into stable beta-N2 crystals.

We report on results on the crystal nucleation and growth of nitrogen. Using molecular dynamics simulations, we show that while nucleation proceeds into the metastable alpha-phase (i.e., the crystalline phase associated with the lowest free energy barrier of formation), growth of the crystallite proceeds through a reorganization of the nucleus into the thermodynamically stable beta-phase.